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UNCLASSIFIED 


Coherent  and  Incoherent  Phase  Retrieval  using 

Gaussian  Beam  Modes 


Creidhe  O'Sullivan1,  Stafford  Withington2  and  J.  Anthony  Murphy1 


Abstract  - Phase  retrieval  techniques  are  important  at 
terahertz  frequencies  where  it  is  often  difficult  to  determine 
the  phase  of  a field  directly.  In  this  paper  we  present  a 
phase  retrieval  technique  based  on  a Gaussian-Beam 
formalism.  We  show  how  the  complex  mode  coefficients  of 
a coherent  field  can  be  determined  by  fitting  to  intensity 
distributions  at  any  two  planes.  We  extend  the  analysis  to 
incoherent  or  partially-coherent  fields  by  working  in  terms 
of  coherence  matrices  rather  than  mode  coefficients. 
Working  with  the  coherence  matrix  also  allows  us  to  vary 
the  entropy  of  the  field  as  part  of  the  fitting  process  and  to 
handle  the  problem  of  noise  in  the  measurements  more 
correctly.  We  illustrate  our  techniques  with  some  simple 
examples.  ’ 


I.  Introduction 

The  phase  retrieval  problem  arises  when  attempting  to 
determine  the  phase  of  a complex  function  from 
measurements  of  intensity  alone.  Such  situations  occur 
frequently  in  fields  such  as  microscopy  and  antenna 
design  where  far-field  intensity  patterns  can  most  easily 
be  measured.  The  complex  coefficients  of  non-ideal 
horns  and  of  antennas  ( e.g . planar  antennas)  whose  mode 
coefficients  are  not  otherwise  known,  are  needed  in  the 
design  of  many  submillimetre  optical  systems.  Different 
object-plane  fields  can  yield  the  same  intensity  in  the  far- 
field,  and  so  to  obtain  a unique  solution  for  the  object- 
plane  phase,  additional  information  must  be  supplied  (see 
eg.  Taylor[l]  for  an  overview  of  the  problem  and 
possible  solutions).  Gerchberg  and  Saxton  [2]  suggested 
a technique  using  intensity  measurements  in  the 
Fraunhofer  and  object  planes,  Misell  [3]  used  intensity 
distributions  in  two  slightly  defocussed  images.  In  the 
case  of  few-moded  bolometers  for  example,  the  analysis 
must  be  applicable  to  partially  coherent  fields. 

Here  we  follow  the  method  of  Isaak  et  al.[4]  and  describe 
the  object  in  terms  of  Gaussian  beam  modes  (see  e.g. 
[5]).  The  mean-sqaure  error  between  the  measured 
(numerically  simulated)  and  trial  intensity  distributions 
on  two  planes  is  then  minimised  by  adjusting  the 
complex  mode  coefficients.  Again  using  the  Gaussian 
beam  formalism,  we  extend  the  technique  to  partially 
coherent  fields  by  adjusting  the  elements  of  its  coherence 
matrix  rather  than  its  mode  coefficients.  We  find  it 
useful  to  describe  the  degree  of  coherence  of  the  field, 
using  the  concept  of  entropy,  in  order  to  constrain  the 
solution,  particularly  for  the  incoherent  case  or  where 
noise  is  present  in  naturally  coherent  fields. 
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II.  Coherent  Phase  Retrieval 


The  intensity  distribution  of  a beam,  at  a distance  z from 
the  waist,  can  be  expressed  in  terms  of  basis  functions 
(y/n)  appropriate  to  the  symmetry  of  the  source  as 


I(r,z)a c 


exp  (jPn  ))vn  </>  *)  exp(/  «AfzSO)) 


, (1) 


where  an  and  pn  are  the  mode  amplitudes  and  phases, 
the  phase  slippage  between  modes  is 

A</>  = tan_1f— :L— 1 , (2) 

\nW  J 

and  W is  the  beam  waist  radius. 


Table  1:  Coherent  phase  retrieval 


Actual  Values 

Retrieved  Values 

a, 

A Pi 

Q\ 

bpi 

0.78119 

- 

0.78125 

- 

0 

(0.1000) 

1.05055x10^ 

-0.6233 

-0.23166 

0.2000 

-0.23177 

0.1997 

0 

(0.3000) 

-4.25362x1  O'4 

3.0074 

0.053578 

0.4000 

0.053841 

0.4092 

So  the  phase  retrieval  problem,  in  the  coherent  case,  is  to 
determine  the  an  and  the  pn  by  fitting  to  intensity 
measurements  alone.  Intensity  measurements  on  at  least 
two  planes  are  needed,  but  unlike  some  retrieval 
techniques  there  is  no  restriction  on  which  two,  since  the 
beam  can  be  propagated  to  any  plane  by  simply  changing 
the  phase  slippage  term.  The  expression 

(3) 

k 

is  minimised  using  a suitable  algorithm1"  to  adjust  the 
mode  amplitude  and  phases.  4,  is  the  modelled  intensity 
at  point  r,-  on  plane  k and  is  the  corresponding 
measurement.  (Sampling  and  data-weighting  issues  have 
been  discussed  elsewhere  [4]).  As  an  example  we  have 
simulated  the  (one-dimensional)  field  produced  by  a 
corrugated  conical  horn  at  two  planes,  one  close  to  the 
mouth  of  the  horn,  the  other  in  the  far-field.  We  sampled 


+ Minimisations  were  carried  out  using  MatHEMATIA®’s 
FindMinimum  routine  which  uses  a modification  of  Powell's  method 
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the  distributions  at  1 1 points  and  attempted  to  recover  the 
first  5 Hermite-polynomial  modes. 

The  actual  values  used  in  the  simulation  are  listed  in 
Table  1 along  with  the  retrieved  values.  The  initial  trial 
values  were  zero,  except  for  the  lowest-order  mode 
amplitude  which  was  set  to  1.  We  found  the  coherent 
phase  retrieval  for  examples  of  this  type  did  not  depend 
significantly  on  the  initial  values  chosen.  The  phase 
retrieval  process  does  not  constrain  the  absolute  phase 
and  so  we  list  the  phase  difference  between  each  mode 
and  the  fundamental. 

Table  1 shows  excellent  agreement  between  retrieved  and 
simulated  values.  Optimisation  of  the  choice  of  planes, 
measurement  weightings  and  sampling  may  further 
improve  the  results,  but  in  this  paper  we  instead 
concentrate  on  extending  the  technique  to  partially 
coherent  fields.  As  the  method  described  so  far  relies  on 
the  field  modes  maintaining  a fixed  phase  relationship, 
another  method  must  be  chosen  for  this. 

An  alternative  approach,  appropriate  for  both  coherent 
and  partially-coherent  fields,  is  to  work  in  terms  of  their 
coherence  matrices  rather  than  mode  vectors.  We 
describe  such  a method  in  the  following  sections. 

III.  The  Coherence  Matrix 

The  analysis  of  partially  coherent  fields  has  been 
described  previously  [6].  We  summarise  some  of  the 
main  points  here.  A partially  coherent  field  is 
constructed  from  a set  of  coherent  diffracting  free-space 
modes.  We  assume  that  the  field  under  investigation  is 
one  member  of  an  ensemble  of  such  fields  and,  as  with  a 
coherent  field,  it  can  be  expanded  in  terms  of  a sum  of 
modes 

E ' (r> ©)  = 2 A‘m  (®V/»  (r> ®) > (4) 

m 

so  long  as  the  bandwidth  is  sufficiently  small  that  the 
phase  at  one  point  in  a member  of  the  ensemble  is  well- 
defined  in  respect  to  the  phase  at  any  other  point  in  the 
same  field.  The  cross  spectral  density  W is  then 

W(r'>r)  = ^(r)Er(r'jj 

= X ('■'M'Om  (5) 

m n 

where 

cmn=(4  4)  (6) 

are  the  elements  of  the  coherence  matrix.  The  coherence 
matrix  characterises  the  form  of  a field  at  any  plane,  with 
all  the  second-order  statistical  properties  being 
completely  specified.  In  terms  of  the  cross  spectral 
density,  the  elements  of  the  coherence  matrix  are  given 
by 

Cmn  = J [W(r\ r)j/n  (r)ds ds'  (7) 

where  the  expansion  functions  y/m  form  an  orthonormal 
basis  set  and  5 is  the  source. 

For  an  incoherent  field  of  uniform  intensity 


W (/■', r)=  I0S(r' - r)  and  (8) 

c = V;  (9) 

all  the  modes  are  excited  equally  and  independently.  For 
a completely  coherent  field  the  coherence  matrix  can  be 
simply  calculated  from  the  ordinary  mode  coefficients  as 

C = AA*T  » (10) 

A partially  coherent  field  can  be  traced  through  a 
submillimeter-wave  optical  system  using  the  overall 

coherent-mode  scattering  matrix  S of  the  system. 
Withington  & Murphy[6]  have  shown  that 

D = SCS*r  (11) 

where  C is  the  coherence  matrix  at  the  input  plane  of  the 
system  and  D is  the  coherence  matrix  at  the  output  plane. 
The  elements  of  the  scattering  matrix  for  free-space 
propagation  to  a plane  a distance  z away  are  simply 

sm„  = expO'mA^K™  (12) 

where  is  again  given  by  (2). 

Just  as  in  the  coherent  case  where  the  fields  were 
decomposed  into  modes  and  then  propagated  from  one 
plane  to  another,  in  the  more  general  case  the  cross 
spectral  density  is  broken  down  into  modes  and  it  is  these 
modes  (or  elements  of  the  coherence  matrix)  that  are 
propagated.  The  intensity  distribution  of  the  field  can 
then  be  recovered  from  the  coherence  matrix  using 

l(r)  = II  c.  • 03) 

m n 

In  phase  retrieval  process  we  determine  those  elements  of 
the  coherence  matrix  that  give  the  best  fit  to  two 
measured  intensity  distributions.  The  NxN  coherence 
matrix  of  a fully  coherent  field  (given  by  (10))  can  be 
constructed  from  2AM  independent  quantities 
corresponding  to  the  N mode  amplitudes  and  AM  phase 
differences  of  §1.  Phase  retrieval  results  using  the 
coherence  matrix  formalism  were  similar  to  the  mode 
vector  method,  although  it  was  found  that  for  some  initial 
trial  values  the  minimisation  routine  did  not  converge. 


(a)  (b) 


Fig.  1:  Intensity  (arbitrary  units)  as  a function  of 

off-axis  distance  for  an  incoherent  beam  in 
(a)  the  near-  (continuous  line)  and  far-field 
(dashed)  and  (b)  in  the  source  plane  (the 
continuous  line  is  the  recovered  intensity).  5 
modes  were  used  and  the  intensity  was 
sampled  at  21  points  along  the  two  planes 
plotted  in  (a). 

For  an  incoherent  or  partially  coherent  field,  however,  the 
only  restriction  we  have  placed  on  the  matrix  is  that  it  is 
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Fig.  2:  Intensity  (arbitrary  units)  as  a function  of  off  axis 

distance  for  a beam  in  the  (a)  near-  and  (b)  far-field 
and  in  the  (c)&(d)  object  plane.  The  continuous  line 
shows  the  intensity  derived  from  the  recovered 
coherence  matrix  and  the  dashed  line  the  actual 
intensity.  5 modes  were  used  and  the  intensity  was 
sampled  at  21  points  along  two  planes. 


Hermitian,  and  the  phase  retrieval  fitting  process  was 
carried  out  with  respect  to  N2  variables. 

The  coherence  matrix  of  a fully  incoherent  source  with  an 
intensity  distribution  I(r)  can  be  calculated  using  (7)  with 
W(r',r ) = /(rj^r'-/1).  (Despite  the  fact  that  the  source 
itself  is  incoherent,  correlations  will  exist  between  modes 
if  the  intensity  distribution  is  not  uniform  [6]).  We 
simulated  an  incoherent  source  with  the  same  intensity 
distribution  as  the  example  in  §1  and  carried  out  the 
retrieval  process  for  an  Hermitian  coherence  matrix.  The 
simulated  intensities  in  the  near  -and  far-field,  along  with 
the  retrieved  and  simulated  intensities  in  the  source  plane 
are  plotted  in  Fig.  1.  The  elements  of  the  retrieved 
coherence  matrix  differed  from  the  actual  values  by  a few 
percent. 

Surprisingly,  the  success  of  this  technique  was  found  to 
decrease  with  increasing  beam  coherence.  Fig.2 
illustrates  the  results  for  a completely  coherent  beam 
(again  making  no  restriction  on  the  coherence  matrix 
other  than  it  be  Hermitian).  A coherence  matrix  that 
described  the  field  well  at  the  two  sampled  planes  was 
found  and  yet  this  solution  was  not  a good  approximation 
to  the  original  coherent  field.  Quite  often,  if  the  solution 
was  not  constrained  to  be  a coherent  field  the  method 
converged  to  an  incoherent  one.  Sampling  at  more  than 
two  planes  improved  the  result  but  the  since  the  fit  to  the 
measured  data  was  always  very  good  there  was  no  way  of 
predicting  the  accuracy  of  the  reconstructed  object  field. 

So,  for  completely  coherent  or  incoherent  fields  we  have 
retrieval  methods  that  recover  good  approximations  to  the 
coherence  matrix.  For  partially  coherent  fields,  however, 
the  minimisation  must  converge  to  the  correct  solution 
which  has  a beam  coherence  somewhere  between  the  two 
extremes. 

For  this  we  need  to  supply  information  that  characterises 
the  coherence  of  the  beam.  A coherent  field  has  a single 
dominant  natural  mode  whereas  for  an  incoherent  field 


power  is  spread  evenly  amongst  each  of  the  natural 
modes  of  the  field.  One  way  of  describing  this 
distribution  of  power  amongst  field  modes  is  by  the 
entropy  of  the  field. 


IV.  Entropy 

It  turns  out  that  the  entropy  of  the  field  can  be  easily 
determined  once  the  coherence  matrix  is  known[9]. 
Since  the  NxN  coherence  matrix  of  the  field,  C,  is  non- 
negative and  Hermitian,  it  has  N real  eigenvalues 
Each  eigenvalue  corresponds  to  the  total  power  in  each 
natural  mode  of  the  field.  The  probability,  ph  that  a 
detected  photon  belongs  to  the  P eigenmode  is  given  by 

Pi  ~ Y 7 ’ CD/*/  = 0 04) 

j 

and  the  entropy  (H)  of  the  system  becomes 


H = ~YJpMpi  (15) 

1=1 

where  the  sum  is  taken  over  all  possible  modes  (TV). 

A completely  coherent  wavefield  is  specified  by  only  one 
non-vanishing  eigenvalue,  whereas  a completely 
incoherent  field  (maximum  amount  of  disorder)  has  all 
eigenvalues  the  same.  The  entropy  of  a field  can  thus 
vary  between  0 (fully  coherent)  and  ln(jV)  (fully 
incoherent).  As  an  example  The  entropy  of  the  field  in 
Fig.  2 was  1.06  (compared  to  a maximum  possible  of 
1.609)  which  illustrates  the  coherence  imposed  on  the 
field  due  to  the  non-uniform  intensity  distribution. 
Diagonalising  the  coherence  matrix  and  finding  its 
eigenvalues  allows  the  entropy  of  the  field  to  be 
determined. 

Repeated  diagonalisation  of  the  coherence  matrix  in  a 
minimisation  routine,  however,  can  be  computationally 
expensive  and  so  we  use  an  alternative  way  of 
determining  the  entropy  from  the  coherence  matrix. 

It  can  be  shown  that  the  entropy  of  a wavefield  can  also 
be  expressed  in  terms  of  its  coherence  matrix  as 


H = -7V[C' ln(C')] 


(16) 


or,  using  the  power  series  expansion  of  natural  log 


(17) 


where  / is  the  identity  matrix  and  Cr  is  the  normalised 
coherence  matrix  (C/Tr[C]). 


The  value  of  the  entropy  obtained  using  the  above  power 
series  converges  reasonably  quickly  to  the  value 
calculated  by  means  of  diagonalising  the  coherence 
matrix.  To  check  the  rate  of  convergence  we  have  used 
both  methods  to  calculate  the  entropy  of  a well-behaved 
source  of  variable  entropy.  A Gauss-Schell  source[7]  has 
a Gaussian  intensity  distribution  (width  characterised  by 
<rs)  and  a Gaussian  coherence  length  (characterised  by 
o*g).  By  changing  the  relative  scale  factors  (< ^ it  is 
possible  to  move  from  having  a fully  coherent  Gaussian 
through  to  having  a fully  incoherent  Gaussian.  We 
calculated  the  entropy  of  the  source  as  a function  of  the 


degree  of  coherence  using  both  the  diagonal isation  and 
power  series  methods. 
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Fig.  3:  Plot  of  entropy  as  a function  of  the 

number  of  terms  in  the  expansion  of 
equation  (17).  Entropy  is  calculated 
for  four  Gauss-Schell  sources  with 
different  relative  scale  factors  as/ag  (ao 
for  a fully  incoherent  source,  0 for  a 
fully  coherent  source).  The  dashed 
lines  show  the  entropy  calculated  by 
diagonalising  the  coherence  matrix. 
(Since  only  5 modes  were  used  in 
these  calculations,  the  maximum 
entropy  is  loge5  = 1.609) 

Fig.  3 shows  that  the  power  series  expansion  is  a good 
way  of  getting  a first-order  measure  of  the  degree  of 
overall  coherence  efficiently. 

V.  Partially  Coherent  Fields 

Finally  we  look  at  an  example  of  a partially  coherent 
field.  We  have  chosen  a fully  coherent  horn  surrounded 
by  a sheet  of  perfect  absorber  (a  fully  incoherent  source). 
The  coherence  matrix  of  the  whole  system  is  simply  the 
sum  of  the  individual  coherence  matrices  for  the  horn  and 
absorber  [6].  Fig.  4 shows  the  result  of  starting  with  a 
uniform,  fully  incoherent  source  and  minimising  both  the 
intensity  difference  and  the  entropy  of  the  field.  Without 
the  use  of  entropy  an  excellent  fit  to  the  data  was  found, 
(total  error  -1  O'9)  but  the  recovered  field  did  not  resemble 
the  actual  source  field.  Minimising  the  entropy  as  well  as 
the  intensity  error  gave  a much  improved  estimate  for  the 
original  coherence  matrix. 

VI.  Conclusion 

This  paper  describes  numerical  simulations  to  test  phase 
retrieval  using  a Gaussian-beam  formalism  to  represent 
complex  fields.  For  coherent  fields  fitting  for  mode 
vectors  was  found  to  be  a robust  technique.  We  have 
tested  both  coherent  and  partially  coherent  fields  using  a 
coherence  matrix  method.  As  with  other  techniques, 
success  depends  on  the  function  being  minimised, 
sampling  and  the  initial  estimate,  but  one  significant 
advantage  with  this  method  is  that  it  does  allow  intensity 
measurements  to  be  taken  across  any  number  of  planes. 
Here  we  have  restricted  ourselves  to  two.  For  coherent 
fields  it  was  found  that  local  minima  did  not  vary  too 
much  from  the  true  value,  but  the  minimisation  did  not 
converge  for  some  initial  trial  values.  For  partially 


distanceffom  axis 

Fig.  4:  Intensity  (arbitrary  units)  as  a function 

of  off-axis  distance  for  a coherent  horn 
surrounded  by  an  incoherent  absorber. 
The  continuous  line  shows  the 
recovered  intensity,  the  broken  line  is 
the  actual  intensity.  5 modes  were 
used  to  model  the  field,  and  the 
intensity  was  sampled  at  21  points 
along  two  planes. 

coherent  fields  local  minima  were  often  found  that  fitted 
the  simulated  data  very  well  but  did  not  reproduce  the 
object  field.  We  found  that  incorporating  the  entropy  of 
the  beam  expected  into  the  function  to  be  minimised 
improved  the  results. 

Future  work  will  look  at  the  effect  of  measurement  error 
on  the  phase  retrieval  technique,  as  well  as  a closer 
investigation  of  sampling,  weighting  of  data  and  the  best 
function  to  minimise. 
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